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Abstract. The difficulty of description of the radiative transfer in disordered 
photonic crystals arises from the necessity to consider on the equal footing the 
wave scattering by periodic modulations of the dielectric function and by its 
random inhomogeneities. We resolve this difficulty by approaching this problem 
from the standpoint of the general multiple scattering theory in media with 
arbitrary regular profile of the dielectric function. We use the general asymptotic 
solution of the Bethe-Salpeter equation in order to show that for a sufficiently 
weak disorder the diffusion limit in disordered photonic crystals is presented 
by incoherent superpositions of the modes of the ideal structure with weights 
inversely proportional to the respective group velocities. The radiative transfer 
and the diffusion equations are derived as a relaxation of long-scale deviations 
from this limiting distribution. In particular, it is shown that in general the 
diffusion is anisotropic unless the crystal has sufficiently rich symmetry, say, the 
square lattice in 2D or the cubic lattice in 3D. In this case, the diffusion is isotropic 
and only in this case the effect of the disorder can be characterized by the single 
mean-free-path depending on frequency. 
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1. Introduction 

The interest in structures with periodic modulations of the dielectric function 
(photonic crystals lj) has been motivated initially by the possibility to modify 
substantially the spontaneous emission in such media. One of the main motives 
of the study, therefore, has been the existence of the complete band gap 2J when 
the propagation of light is completely inhibited inside some frequency region. Only 
relatively recently it has been realized that the periodicity of refractive index in 
photonic crystals by itself results in a number of unusual properties even in absence 
of the complete band gap [3J |U O [7J [5] ■ These properties do not necessarily require 
a strong contrast of the periodic modulation, i.e. the ratio of the minimum and the 
maximum values of the refractive index, and, particularly, can be observed even if 
the contrast is much weaker than needed for the gap to open. One of the illustrative 
examples is provided by dark modes [9l I10[ lllj , which are not coupled to plane waves 
propagating outside the structure. The existence of these modes is related to the 
point symmetry of the photonic crystal and, therefore, the dark modes present even 
if the contrast is very low. The dark modes, as well as other phenomena such as self- 
collimation and the negative refraction demonstrate that periodic spatial modulations 
of the dielectric function have more ways to affect the propagation of light than just 
producing a gap in the spectrum. 

Real photonic crystal structures always contain one or another type of disorder 
regardless of a manufacturing procedure. It is crucially important, therefore, to 
understand to which extent disorder affects properties of these structures. This issue 
is of a great interest because an interplay between periodic and random variations 
of the refractive index creates new challenges for a theory of light propagation 
in inhomogeneous media, and promises new and unusual effects in the radiative 
transport. 

The problem of the disorder in photonic crystals can be approached from two 
perspectives. On one hand, there is an issue of effects of disorder on spectral features 
of photonic crystals and their manifestation in such characteristics as reflection or 
transmission spectra. This research direction involves, for instance, studying such 
problems as a dependence of the width of the photonic band gap on the degree of the 
disorder [121 E2 03] • A different type of questions arise when one is concerned with 
effects of disorder on propagation of light inside the photonic structure within the 
framework of the transport theory. Problems considered in this case include diffusion 
in the photonic crystals [HI [16] or enhanced backscattering [UJ [TBI Qjj] . 

The main objective of the current paper is to develop a general theoretical 
approach to wave transport in disordered photonic crystals, which would 
systematically, from first principles, take into account the periodic nature of the 
average refractive index. The microscopical nature of our approach distinguishes it 
from earlier papers, which relied either on ad hoc modifications of results obtained 
within the multiple scattering theory in statistically homogeneous media |19| 118) . 
or on phenomenological assumptions regarding the distribution of the field in the 
bulk of the photonic crystal [TBI [TjJ. For instance, it was assumed in the latter 
papers that propagation of light in the bulk of the photonic structure is the same as 
in statistically homogeneous random media, and the band structure of the photonic 
crystal only manifests itself within a narrow surface layer of the sample. This idea is 
obviously based on the assumption that multiple scattering destroys photonic modes 
of an ideal periodic structure so that the structure of the electromagnetic field in 
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the bulk of the photonic crystal becomes indistinguishable from that of a regular 
disordered medium. A number of experimental and numerical results, however, cast 
doubts on this assumption. For instance, it was shown in [20) that photonic modes are 
completely developed in samples with linear dimensions as small as just a few periods. 
It was also determined experimentally that the mean- free-path, £, [see Eqs. ((60]) and 
l|72p ] due to disorder in real photonic crystals (see. e.g. Table I in [21]) substantially 
exceeds the lattice constant of the photonic crystal. Even for relatively high frequencies 
Vta/2-Kc — 1.6 it was found that IJa ~ 4. Thus, even in this least favourable case there 
are about hundred elementary cells in a volume with the linear dimensions of the order 
of I. Comparing these two results and invoking ideas of the separation of length scales 
it is reasonable to expect that the underlying periodicity of the photonic crystals must 
still manifest itself even in the presence of disorder. In our paper we show that this 
expectation is indeed justified, and demonstrate explicitly the effects of periodicity on 
wave transport in the diffusive regime. 

The standard multiple-scattering theory of wave transport in disordered media 
depends significantly on the plane-wave representation of the scattered field. The 
role of the plane waves is explicitly emphasized by the use of the Wigner function 
in derivations of the radiative energy transfer equations [22 , 23l [24l [25] . In photonic 
crystals, however, approaches based on plane waves encounter significant difficulties 
because plane waves are not normal modes of the underlying ideal periodic structure. 
While the Wigner representation of the field-field correlation function itself remains, 
of course, valid even in such structures, the Wigner function, however, ceases to be a 
smooth function of coordinates, which is essential for the derivation of the radiative 
transfer and diffusion equations. 

Qualitatively, one of the difficulties of the plane wave based approaches is due 
to the fact that the plane waves are scattered not only by the random fluctuations 
of the refractive index but also by its periodic modulation. The latter is a purely 
deterministic process and is responsible for the formation of photonic crystal modes, 
which can be considered as coherent superpositions of the plane waves. Thus, in 
order to describe wave transport in disordered photonic crystals one has to be able 
to separate the deterministic contribution to the scattering from the one caused 
by the disorder. This can be achieved by developing the transport theory on the 
base of photonic modes of an ideal crystal, but as we will see even within this 
approach the discrimination between coherent and incoherent processes is highly 
non-trivial. The second, more technical, difficulty in adapting the standard multiple 
scattering theory to photonic crystals arises from the fact that the theory developed 
for statistically uniform media heavily depends on certain assumptions (e.g. the 
translational invariance of the averaged Green's function) that are not valid in 
disordered photonic crystals. 

In the present paper we resolve these difficulties and develop a consistent multiple 
scattering theory of light transport in media with the periodic-on-average dielectric 
function. Some of the results obtained are actually valid also in media with arbitrary 
modulation of the background (average) dielectric functions. We introduce the 
interpretation of the field-field correlation function as the density matrix and show 
how it can be used to separate coherent and incoherent contributions in the transport. 
Using this idea we generalize the concept of specific intensities to the case of photonic 
crystals and derive respective radiative transfer equations. We also find an asymptotic 
solution of Bethe-Salpeter equation describing a steady state intensity distribution in 
an infinite media far away from sources, which we use to derive a diffusion equation in 
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a steady state regime describing long scale spatial relaxation of the intensity toward 
the limiting distribution. 

2. Multiple scattering in disordered photonic crystals 

In the framework of the scalar model, the spatial distribution of the wave field at 
frequency u in a disordered photonic crystal is governed by the Helmholtz equation 

AE u (r)+u 2 e(r)E ul (r)=j(r), (1) 

where j(r) is the external source. Here and henceforth we use units with c = 1. In (fTJ) 
we have introduced the dielectric function 

e(r) = e(r) + Ae(r), (2) 

which consists of two components. The periodic part, e(r + a) = e(r), constitutes the 
photonic crystal with a being the vector of lattice translations. The zero-mean random 
term, Ae(r), describes the deviation of the dielectric function from the ideal periodic 
form. We assume that the random part of the refractive index is a zero-mean Gaussian 
random field, i.e. its statistical properties are completely characterized by the 
covariance K{r\,r2) = (Ae(ri)Ae(r2)}. In |Appcndix A| we provide a model of such 
inhomogeneities relevant for disordered photonic crystals. With this assumption (fTJ) 
can be readily analyzed by standard diagrammatic techniques based on the Born 
series representation of the solutions of the integral (Lippmann-Schwinger) formulation 
of ([T]) , which were developed in the theory of multiple wave scattering in statistically 
uniform media [26[ I27[ 128] . This approach is virtually model independent and can be 
applied for structures with arbitrary spatial profile of the deterministic part of the 
refractive index. 

The Dyson equation for the Green's function of JT]) averaged over realizations of 
the disorder, G = (G), is obtained using the standard diagrammatic technique and 
has the form 

G ul (r,r') = G (r,r') + J dndr 2 G (r, njE^n, r 2 )G u {r 2 , r'), (3) 

where E w is the self-energy "defined" as a sum of all irreducible diagrams. We use 
somewhat cumbersome coordinate representation in order to emphasize its generality 
and independence of a particular form of the regular modulation e(r). The latter 
determines the nonperturbed Green's function, Go(r,r'), which is assumed to be 
known. In particular, in disordered photonic crystals Go(r,r') is the Green's function 
of the ideal periodic structure. 

The main difference between statistically uniform media and disordered photonic 
crystals is reflected in this equation through symmetry properties of the unperturbed 
Green's function. In periodic-on-average systems, this function is invariant with 
respect to lattice translations and the group of point symmetries of the underlying 
photonic structure, while in the uniform disordered systems it has full translational 
and rotational symmetry. In | Appendix B| we demonstrate that the averaged Green's 
function and the self-energy posses the same translational and point symmetries as 
Go(r,r'). We utilize these properties by expanding all related quantities in terms of 
photonic Bloch modes, 

*„,fc(r) = j k - r u k , n (r), (4) 

where u n> k{r) is periodic with the period of the photonic lattice, fe is the Bloch wave 
vector lying inside the first Brillouin zone, and n enumerates photonic bands. Using 
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these functions we introduce the matrix representations of the quantities appearing 
in © according to 

G (n,r 2 ) = ^Sn( fe )*fc,n(ri)tf; in (r 2 ), 

h . n 

G w (r ls ra)= ^ ^W^mK)*^^), 

h,n,m 

E w (n,r 2 ) = ^ E nim (fc)e(ri)tf fc , n (ri)*i ifn (r 2 )e(r2). 
The respective matrix elements in these expansions are 



(5) 



W 2 -£j2(fc)' 

G n , m (fe) = J dr i^ r 2e{ri)u* k n (r 1 )G^(r 1 ,r2)u k ^ m (r 2 )e(r 2 ), (6) 

E n , m (fc) = ^5 J ^ r i Ar 2U kn {ri)Y, u {r 1 ,r2)u ktm {r 2 ), 

where w^(fc) is the dispersion law of n-th band, the integration is performed over the 
elementary cell of the ideal structure, and V is the volume of the elementary cell. 
Deriving ([6]) we use the orthogonality condition of the Bloch functions 

J dre(r)*;, fc (r)* m ,,(r) = S mn 6(k - q). (7) 

Using ^ we can rewrite ([3]) in the matrix form 

Gmn(k) = g m (k)S mn + y^g m (fc)S m ;(fc)g; n (fc). (8) 

i 

This equation emphasizes the fact that the translational invariance of the self-energy 
prevents modes with different Bloch vectors to be mixed while states corresponding to 
the same Bloch vector but belonging to different bands are coupled by non-diagonal 
elements of the self-energy E m j(fe). 

In order to analyze the general effect of the disorder, we separate the diagonal, 
Ed(fc), and the off-diagonal, E (fe), parts of the self-energy representing the latter in 
the form 

E(fe) = E d (fc) + E (fc). (9) 

The diagonal part Ed(fc) modifies each band independently. It can be accounted for 
by introducing a modified Green's function 

Go(fc) = : ^ , (10) 

Go X (fc)-Ed(fc) 

which is determined by the Dyson equation similar to the one written for the standard 
case of a statistically homogeneous medium. In terms of the modified Green's 
function © takes the form 

G = G + G Z o G. (11) 

Eqs. pop and (fTTj) show the two-fold role of the disorder in disordered photonic 
crystals. The disorder not only modifies each band separately similar to the case 
of the statistically homogeneous media, but also couples these modified bands. It 
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is important to note that the band coupling is a subject of various selection rules. 
First, the most restrictive rule comes from the translational symmetry of the self- 
energy. As has been noted, it prevents states characterized by different Bloch vectors 
to be coupled. In other words, from the perspective of a band diagram one can have 
only "vertical" coupling. The second rule follows from the point symmetries. The 
self-energy transforms according to identity representation of the symmetry group 
of a given point in the reciprocal space of the photonic crystal. As a result, its 
matrix elements between states corresponding to different irreducible representations 
vanish. This selection rule has important implications for high-symmetry points and 
directions. In particular, it means that disorder does not lift the degeneracy at the 
points where the degenerate states are described by different or by multidimensional 
irreducible representations. The first can be shown by direct calculations of the 
respective matrix elements of the self-energy. The second follows from the following 
argument. The states corresponding to a representation with a dimension higher 
than 1 can be coupled only with the states that transform according to the same 
presentation. As a result the modified state also transforms according to this 
presentation and, hence, the respective state should also be degenerate. This implies, 
in particular, that the degeneracy is not lifted at high-frequency T-points. 

In order to analyze the effect of the band coupling in more details, we consider a 
two-band model when all bands but two (denoted by 1 and 2) remain uncoupled. The 
solution of the Dyson equation describing the coupled bands has the form 

= W~ ( * 'e" S " a-^\ )' < 12 > 

L>\1 \ ^21 9i ~ Mi J 

where gi and Ey are the respective matrix elements of the unperturbed Green's 
function and the self-energy, respectively. The zeros of the function 

£>i2 - (gi 1 - En) (g^ 1 - £ 22 ) - Si 2 E 2 i (13) 

give the spectrum of averaged excitations. Introducing uji. 2 (k), the unperturbed 
dispersion laws of the interacting bands, the poles of averaged Green's function can 
be written in the form 

^ = I(^+s 2 2 )±^v/rT^. (w) 

Here u)\ 2 = Wi 2 ~ ^11,22 are the band frequencies after the diagonal modification, 
A12 =u)\ — a>2, and 

Ei 2 E 2 i . . 

*?12 = 4 A2 ( 15 ) 
^12 

is the parameter characterizing the strength of the band coupling. As one could expect, 
this parameter is proportional to the matrix elements of the self-energy between 
the coupled bands and is inversely proportional to the frequency separation of the 
modified bands. When this parameter is small, 7712 <C 1, the effect of the band 
coupling can be neglected and the photonic crystals can be described in a single band 
approximation. It should be noted, however, that in spite of formal similarity between 
the averaged Green's function in the single band approximation and the respective 
expression for the uniform disordered medium, the transport properties of the two 
systems remain substantially different. The condition 7712 <C 1 can be easier satisfied 
at lower frequencies, when the band separation is of the order of magnitude of the 
fundamental band-gap. For higher frequency bands, however, the band coupling may 
play a significant role even in the case of photonic crystals with weak disorder and 
strong contrast of the refractive index. 
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3. Transport in disordered photonic crystals 

3.1. General formalism 

The transport properties of disordered media are characterized by the field-field 
correlation function 

AWi,r 2 ) = {E^{ ri )E^{r 2 )) , (16) 

which can be used to describe transfer of energy and, generally, spatial distributions 
and time evolutions of any quantity quadratic in the field. A relation between 
Pwi.w 2 ( r ii r 2) an( i the external sources is provided by the intensity propagator 
U{r 1 ,r 2 ;r' 1 r' 2 ) = (G u)1 (r x ,r' x )Gl 2 {r 2 ,r' 2 )) according to 

Pu 1 ,w 2 (ri,r 2 ) = J dr' 1 dr 2 U(r 1 ,r' 1 ;r 2 ,r 2 )j L , 1 (r' 1 )j* 2 (r 2 ). (17) 

Using the standard diagrammatic technique one can show that the intensity 
propagator satisfies the Bethe-Salpeter equation, which can be written as 
n(n, r 2 ; r' x r' 2 ) = G m {r x ,r[)G* U2 (r 2 ,r' 2 )+ 

J dr 3 dr' 3 dr 4 dr' 4 G Ui (r x , r 3 )G* W2 [r 2 ,r A )U Uu0J2 (r 3 , r 4 ; r 3 , r' A ) (18) 

xU(r 3 ,r' S r' x r' 2 ). 

The kernel U UllU3 is the irreducible vertex presented formally as a sum of 
irreducible diagrams [26J. We would like to emphasize that similarly to the 
Dyson equation, the Bethe-Salpeter equation holds for an arbitrary regular spatial 
modulation of the dielectric function, not necessarily periodic. We show in |Appcndix B| 
that regardless of the spatial dependence of the regular part of the dielectric function, 
the irreducible vertex possesses an important property of reciprocity 

U UJl . UJ2 {r 1 ,r 2l r' 11 r' 2 ) = U Uuul3 (r' x ,r' 2 ;rx,r 2 ). (19) 

Additionally, in disordered photonic crystals this quantity is invariant with respect to 
lattice translations (see | Appendix Bp suggesting the following representation for the 
vortex 

U uu ^[ri,r2\r' x y 2 ) = e ( r i) e ( r 2)eW)e(r 2 ) 

i>i,2,Zi,2,q; 

Uhfr ill, m\ <?3, <?4)<5 (<?1 + 94 - <?2 - 93) 
*ft,* 1 (^)*«,I.W)**» w (' , j)*;„l 1 (K). 

where the bar over a vector denotes the vector reduced to the first Brillouin zone. 

Using Bethe-Salpeter equation, Eqs. (fT%|) and (TIT)) one can derive an equation 
for the field-field correlation function jO tl)1 , W2 (ri, r 2 ). We present it in an integro- 
diffcrcntial form, which is the most convenient for further analysis. To derive such an 
equation one can apply operators 

^) Al ' 2 + ^ 2 ' (21) 
where indices 1 and 2 indicate a coordinate acted upon by the Laplacian, to both sides 
of (fT7|) . As a result, in the region of space free of external sources, one has 

1 A 1 A 22 

-r^ Al - -^A 2 + w\ - uj a 2 
e(ri) e(r 2 



Puj!.Lj 2 (ri,r 2 ) = 
f (22) 
/ dr[dr' 2 T Ul itl)a (r x ,r 2 ;r' x , r' 2 ) Pull tU>2 (r[ ,r' 2 ), 
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where 

Fmwirur^r'^r'z) = -^-yS wl (n, ri)^(r 2 -r^,) - ^-y£* 2 (r 2 , r 2 )£(ri - r[ 



(23) 



e(n) ^ " e(r 2 ) 

x U UJl ^ (JJ2 (r 1 , r 2 ; r 1 , r 2 ). 

This equation has been a subject of numerous investigations [26l [29], most of 
which were concerned with the radiative transfer or diffusion regimes in statistically 
homogeneous media. In many of those works, diffusion was understood as a transport 
process characterized by asymptotically slow (both in time and space) changes of the 
field intensity. In the spectral domain this behaviour manifests itself in the form of 
the characteristic "diffusion" pole of the intensity propagator, which is proportional to 
(iil — -DQ 2 ) in the limits — > 0, Q — > 0. Here frequency and wave vector transfers, 
Q = lui — ll>2 and Q = q± — q 2 , characterize slow spatiotemporal dynamics of intensity, 
and wave vectors q± and q 2 arise from the plane-wave representation of the scattered 
field. This reliance on the plane waves significantly complicates generalization of 
standard microscopic derivations of radiative transfer or diffusive equations to the 
case of disordered photonic crystals. 

In order to better recognize the source of these difficulties and find a way to 
circumvent them, it is necessary to re-examine basic physical ideas about diffusion 
of light in disordered media. As a first step in this direction, in the current paper 
we consider time independent spatial distribution of the wave intensity in an infinite 
medium far away from the sources. This regime arises in the case of a monochromatic 
source, when the field in the structure harmonically depends on time, oc exp(— iu>t), 
so that f2 = 0. In this case, the equation for the field-field correlation function 
p = {E^E*} (hereafter we omit the lower index corresponding to frequency) can 
be obtained from (|22|) by setting w\ = lo 2 = oj; 

-Ai 1 



p(n,r 2 )= /dridri^n.rajri.riXri.ri), (24) 
e(n) e(r 2 J J J 

where T = Tu,u- Using (|24|) . the optical theorem (the Ward identity) can be derived 
[2"5] in the form F{r, r; r[, r' 2 ) = 0, which is especially useful from the technical point 
of view. Integrating this equation over r we obtain the Ward identity in the form 

E(r a> n) -£*(!•!, ra) = / dri'dr 2 ' [G(r' 2 \r") - G*(r",r 2 ')] U(r'{, r 2 '; n, r 2 ). (25) 

The physical picture of the transport in disordered media is developed using 
close relation of the function p(r\, r 2 ) to both transport characteristics and coherence 
properties of the field. The description of transport (e. g. the energy density and the 
flux) is based on the property that the averaged values of any quantity quadratic in 
field can be expressed in terms of convolution of respective operators with p(ri,r 2 ) 
[see Eqs. (|33|) and (j39j) below]. The coherence properties are described considering 
p( T *i, r 2) as the coherence function [30l 131) . 

A quantity, which simultaneously describes such transport related characteristics 
as energy density and flux and coherence properties of the system, is well known in 
quantum statistics. It is called the density matrix [32] . Indeed, on the one hand, 
the density matrix can be used to calculate current and energy densities, and on the 
other hand, the density matrix describes mixed states, which can be characterized 
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as incoherent superpositions of pure or coherent states. Using the density matrix 
analogy one can think of the field-field correlation function as a characteristic of such 
a mixed state of the wave field in disordered media. Separation of coherent and 
incoherent properties of the field would then involve finding states whose incoherent 
superposition would reproduce the field at the level of function p(ri,r 2 ). However, 
before developing this idea any further, we need to demonstrate that the field-field 
correlation function has indeed all the formal properties of the density matrix. 

We expand p(r\, r 2 ) in terms of the eigenstates of the nonperturbed system, say, 
the modes of the ideal photonic crystal, writing 

P(ri, r 2 ) = ]T p^ tf * M (ri)*:(r 2 ), (26) 

where summation over indices p and v enumerating the eigenstates can involve 
integration, ft is easy to see from the definition of the correlation function that 
coefficients p^ %v constitute a Hermitian matrix. This matrix can be diagonalizcd, 
which corresponds to the spectral representation of the statistical operator in quantum 
mechanics, by means of a unitary transformation to another basis 

= B h ^^ (27) 

so that one has (Mercer's theorem [30 | 133 ] ) 

p{ri,r 2 )=Yn^ R {r 1 )%{r 2 ) (28) 



with p R > 0, which follows from the fact that p(ri,r 2 ) is non-negatively defined 
|34j . The diagonal representation of the correlation function can be interpreted as an 
incoherent superposition or mixture of pure or coherent states ^g, which, in turn, as 
follows from (|27|) are coherent superpositions of the states "I^. In order to clarify the 
exact meaning of this expression let us show that the correlation function can be used 
to calculate the energy density of the field and its Poynting vector in much the same 
way as the density matrix is used in quantum statistics. 

The energy density of the field in a steady state can be presented as 

w(r) = ±[iV 2 e(r)\E(r)\ 2 + \VE(r)\ 2 ], (29) 

where e(r) is the total dielectric function including both regular and the random 
components. By introducing the operator 

w = -(V!-V 2 ) 2 /4, (30) 

one can show that the averaged energy density can be expressed in terms of the 
function 

w(ri,r 2 ) = wp(rx,r2) (31) 



(w(R)} =w(R,R), (32) 



as 



where we used the Helmholtz equation to arrive at (|32|) . One can see that this equation 
can be presented in the typical for quantum statistics form as 

(w(R)) = Tr [pw R ] , (33) 

where wr = 5{r\ — R)5(r 2 — R)w. 
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As an example let us consider the case when the density matrix has the form of 
an incoherent superposition of the eigenstates ^ M of the ideal system. Then from ([55]) 
we find 

(w(R))=J2p»MR)i (34) 

where we have introduced the energy density of the /z-th mode 

w,(R) = €, ^^ ri )^;(r 2 )\ ri=r2=R . (35) 

This expression provides a clear explanation of the notion of incoherent superposition. 
Indeed, the average energy of the field in this expression is presented as a sum of 
energies of individual modes ^^(r) as it is expected for the addition of incoherent 
fields as oppose to the sum of the field amplitudes expected for the coherent fields. 

Similar expressions can be obtained for the average value of the Poynting vector, 
S = hi [E*VE — EVE*] /2, which can be calculated using the operator 



which gives 



where 



S = - y (V 1 -V 2 ), (36) 



(S(R)) = S(R,R), (37) 



5(ri,r a ) = 5p(r 1 ,r a ). (38) 

In the case when the density matrix is diagonal in the basis of the eigenfunctions of 
the ideal system we see again that the average Poynting vector is a sum of Poynting 
vectors of each mode, which indicates the absence of any interference effects in the 
superposition of modes ^(r): 

E 



S(R)) = Tr pS R => P ,S,(R), (39) 



where Sr = 5(r\ — i?)5(r 2 — R)S and S^(R) is the distribution of the Poynting vector 
in the p-th. mode 

S p (fl) = S y^ri)y;(r 2 )\ ri=r2=R . (40) 

It is seen from Eqs. (|3"4"|) and (f3"9")) that p^ have the meaning of the weights of the 
incoherent superposition. More generally, if we normalize p(r\, r 2 ) in such a way that 
Tr[p] = 1, we can interpret p s , the eigenvalues of the matrix /0 Mi „, as the distribution 
function in the space of the states "J^. The values of any quantity, quadratic in field 
averaged over the disorder realizations, therefore, can be calculated as the average over 
the distribution function p R in the similar way as it is done in quantum mechanics. 
We would like to emphasize here that the emergence of the incoherent superposition 
in the problem of wave propagation is directly related to averaging over realizations 
of disorder. Without averaging the density matrix defined in (f^5|) would have a form 
p^ v oc a^a u , and matrices of such form have a single non-zero eigenvalue. As a result, 
the sums in (1341) and (|39[) would consist of just one term indicating that p represents 
a pure or coherent state. 

The notion of incoherent superposition expressed by Eqs. (l34|) and (|39|) also allows 
one to provide a physical meaning for the states VPs diagonalizing the density matrix. 
The spatial field distribution in a random media can be in principle presented as 
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a linear combination of functions from any full system: plane waves, Bloch waves, 
etc. The concept of normal modes as well denned spatial distributions of the fields 
that can be excited separately one from another is not very useful here. Indeed the 
distribution of the field in a random media is so complex that at any given frequency it 
is impossible to excite a single mode out of infinitely many degenerate modes. In this 
situation, functions "Jr play a special role as such distributions of the field whose linear 
combination is purely incoherent in the sense of Eqs. ([34| and (J39J - The form of these 
functions is determined by remaining coherence effects in the scattered waves, and 
thus, by using the density matrix formalism we achieve a separation between coherent 
and incoherent contributions to the energy transport in disordered systems. In other 
words, one can say that the functions diagonalizing the density matrix describe the 
modes, which provide the energy transfer. 

In order to illustrate these general idea, let us consider a case of wave scattering in 
a homogeneous random medium. In the case, of an infinite medium and asymptotically 
far away from the sources the field-field correlation function restores its translational 
invariance: p(ri,r 2 ) — > p(ri — r 2 ). Using plane waves as a basis we can rewrite (|26| 
for this particular case as 

p{ri,r 2 ) = J dqdkp un (k)S(k-q). (41) 

One can see that the density matrix in this case is diagonal in the basis of the plane 
waves so that they are responsible for the incoherent transport. 



3.2. The field-field correlation in an infinite photonic crystal: asymptotic behaviour 

In this subsection we consider a solution of the steady-state Bethe-Salpeter equation 
(|24p for the correlation function p{r\,r2) valid in an infinite photonic crystal 
asymptotically far away from the sources. Besides providing an important non-trivial 
example of the application of our formalism, this solution shows how the periodicity 
of the underlying photonic structure affects asymptote of the intensity distribution 
in the disordered structure and provides a starting point for deriving the diffusion 
equation. Using Eqs. (fT9f and ([25)1 one can check that ((24)) is solved by 

P (oo) (n,r 2 ) = [G*(r 2 ,n) -G(r l5 r 2 )] , (42) 

where N is a normalization constant independent of coordinates. Using the reciprocity 
theorem [26] this is rewritten as 

P (oc) (ri,r 2 ) = -Im[G(n,r 2 )] /N. (43) 

This solution is valid in a general case regardless of the specific form of the regular 
modulation of the dielectric function and the distribution of the disorder. In the case 
of statistically uniform media it is reduced to a form found in Refs. [35j and |36j . 
Thus we can consider this function as the asymptotic distribution of the field-field 
correlator, which can be called for this reason an equilibrium density matrix. 

In the case of disordered photonic crystals, function p(°°)(ri, r 2 ) can be expanded 
in terms of normal modes of the underlying periodic structure. According to ((5|) we 
can present this expansion in the following form 

P (oo) (n,r 2 )= ]T P%l(k)* k , n (ri)n, m (r 2 ), (44) 

fc,n,m 
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where 

p%l(k) = [G; hm (k) - G m , n (k)] /2iN 

1 A,,,,., „ ^ (45) 



/ dn dr 2 e(ri)ttfc n (ri)ImG w (ri,r2)ufe, TO (r2)e(r2). 



iVV 2 

This expression shows that in the basis of normal modes of an ideal periodic structure, 
the density matrix is diagonal with respect to the quasi- wave vector k, but is not 
diagonal with respect to the band indices. Following Eqs. (l26|) . ([27]) . and f28]) we 
diagonalized the matrix pm°n{k) by a unitary transformation 

*fe,m = B rfm{k)^ k,n- (46) 

n 

Since the diagonalization procedure involves only band indices and leaves the Bloch 
vector intact, functions ^>k,m{f) can also be presented in the Bloch form similar to ((4]) 

(r) = e ikr u k ^{r). (47) 

Transformation preserves the scalar product because of unitarity so that 

J dr e(r)*; )VB (r)* g , B (r) - 5(fc - g)<W (48) 

Using this property we normalize defining 

N = —tt f dr e(r)Im[G(r,r)], (49) 

where the integration is performed over the elementary cell of the ideal structure. 
Equations (|42|) and (149|) define /c^ 00 ** as the asymptotic form of density matrix, which 
according to Eqs. (|4"1)) and pi?)) is an incoherent superposition of the states ^fe,™ with 
respective weights. 

As was discussed, the eigenvalues of define the probability distribution 

function on the space of the states ^>k,m{ r ) parametrized by the number of band m 
and by the point in the first Brillouin zone k. As follows from the Dyson equation the 
singularities of the averaged Green's function and, respectively, of the eigenvalues of 
pm^n are determined by the dispersion law of the average excitations [see e.g. fI5]l]. For 
more detailed analysis we consider the situation when we can neglect the band coupling 
so that the averaged Green's function is given by §TU\i . In this case, the eigenvalues of 
the density matrix as functions of the quasi- wave vector have the largest values, when k 
obeys the dispersion equation for a given frequency. In other words, these eigenvalues 
reach maximum values on equifrequency surfaces, F m (uj), corresponding to different 
bands of the ideal photonic crystal. The width of these maxima is proportional to 
Im[Sd m (fe)] [ see ©]• If the disorder is weak in the sense of the Ioffe-Regel criterion 
(see below), then for frequencies, which do not lie at band edges, we can obtain 

P (oo) (n,r 2 ) / dk T^m e-I^W'^-^)!^^)^^^), (50) 

2Nuj ^JF m {u) \v m [k)\ 

where the integration runs along the equifrequency surfaces, v m (k) — VfcW m (fc) is the 
group velocity, and 

= J^k)] V (k) = 1 1(fcRi(fc) . (51) 
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Figure 1. (Colour online) The dependence of l/v m (k) along the equifrequency 
surface corresponding to uja/2ir = 0.41. The horizontal plane is the first Brillouin 
zone. Here and below the calculations are presented for a square lattice 2D PC 
made of dielectric cylinders with the contrast of the refractive indices n = 2 and 
with the radius-period ratio r/a = 0.4. 



Here v m (k) is the unit vector along the direction of the group velocity and we have 
expressed the imaginary part of the self-energy in terms of the respective mean-free- 
path £ m (k), which is defined later in (|60[) . 

As follows from ([50]) in the limit of vanishing disorder, which corresponds to 
the on-shell approximation J26[ [29] [37j [38l [39] in the standard theory of transport in 
statistically homogeneous media, the density matrix takes the universal limit 

^(ri.ra) = W dfc ^ * fc , m ( ri )*^ . m (r 2 ). (52) 

The magnitude of the group velocity is, in general, not constant along the 
equifrequency surfaces. As a result, different states are not equally presented in the 
equilibrium distribution, as is seen from Eqs. (|50p and (|52p . In particular, if there are 
flat bands [ID] with low group velocity near the frequency u, then these bands would 
give the main contribution to ■ Another example of a highly inhomogeneous 
distribution is provided by the frequencies when an equifrequency surface touches the 
boundary of the Brillouin zone. In this case, the magnitude of the group velocity 
becomes very low at the points of contact resulting in an increased weight of the 
respective states in the equilibrium distribution (see Figure [1]). 

It should be noted that the on-shell approximation is poorly suited for studying 
the correlation properties of the field. Indeed, the function Pq°°^ (7*1, r 2 ) does not tend 
to any limit as \r\ — r 2 — > 00 while p^°°^ (tt, r 2 ) in (jST)]) obviously vanishes in this 
limit. Thus, for this purpose one has to use (fi2"]) or in the case of a weak disorder the 
simpler version ([50]) . as it has been practically done in [3S], where additionally the 
approximation of spherical equifrequency surfaces was used. 

However, for studying transport properties the on-shell approximation may by 
appropriate because, as follows from Eqs. ([3"2"]) and {37]), the energy density and 
the Poynting vector are determined by the behaviour of p(ri,r 2 ) near the diagonal 
ri = r 2 . We estimate the effect of the exponential term in (|50p by calculating the 
average energy density. The density matrix in the on-shell approximation yields the 
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energy density as a sum of contributions of different modes [see (I34[) ] 



(53) 



where Wk m (R) is the energy density of the mode ^k,m of the ideal photonic crystal. 
The exponential term in (|50p can be shown to modify each contribution by the factor 



In order to evaluate the correction, we note that ujl m (k) is a parameter of the Ioffc- 
Regel type |27] , which is expected to be much larger than unity far from the localization 
regime, so that the deviation of the equilibrium density matrix from (|52[) can be 
neglected in the analysis of transport quantities. 

The problem of special interest is what happens in an immediate vicinity of the 
complete band gap. Formal application of (foT))) gives vanishing mean-free-path for 
the modes lying at the edge of the gap. This means that approximation ([50)) for the 
density matrix is no longer valid and that the modes of the ideal photonic crystal are 
not a good approximation for the modes Such nonperturbative reconstruction of 
the modes constitutes a special problem and is not considered in the present paper. 
Therefore, in what follows we restrict to the frequencies which are not too close to the 
complete gap so that the condition ut m (k) > 1 can always be fulfilled if the magnitude 
of the inhomogeneities is not too high. 

Let us conclude the analysis of the equilibrium distribution by outlining two 
significant differences between the transport in statistically homogeneous media and 
disordered photonic crystals. The most essential difference is that the radiative 
transport in photonic crystals is provided not by plane waves but by the modes, which 
reflect properties of the ideal structure. Indeed, these are the modes of the photonic 
crystals that constitute /^°°)(ri, 1*2) rather than pure plane waves. As a result, even in 
the asymptotic distribution of the energy density the disorder does not wash out the 
features specific for the ideal structure. In the case of weak disorder when the energy 
density can be approximated by (|53|) . one can see explicitly that frequencies near the 
fundamental gap and higher, the spatial profile of Wk, m {R) is highly inhomogeneous 
on the length scale of the elementary cell |41| . This inhomogeneity is stable, to 
some extent, with respect to averaging over the equilibrium distribution function. 
Most clearly this effect should be seen for frequencies near the lower edge of the 
fundamental gap. Indeed, as follows from the variational principle [1] at such 
frequencies the field distribution takes higher values in the regions with higher values of 
the refractive index for each point on the equifrequency surface. The same obviously 
holds for the energy density of each mode. Thus, we can conclude that even the 
long scale asymptotic behaviour of intensity in disordered photonic crystals is highly 
inhomogeneous. This property is unexpected from the point of view of the transport 
in statistically homogeneous media, where the distinctive feature of this limit is the 
homogeneous distribution of the energy density [42] 




The second important feature of the field distribution in photonic crystals is its 
significant anisotropy as a function of direction of the Bloch vector. Anisotropy itself 
is not, of course, the specific feature of disordered photonic crystals. The similar effect 
can be expected in regular anisotropic media, where the group velocity depends on the 
direction of the wave vector. In disordered photonic crystals, however, this anisotropy 
can be very pronounced, especially near the frequencies corresponding to the edges of 
(partial) band gaps. 



1 



(54) 



1 - 
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3.3. The radiative transfer equation 

Applying ([3"Tj) to the equilibrium state described in the previous subsection, 

S(°°\r 1 ,r 2 ) = Sp(°°Hr 1 ,r 2 ), (55) 

and using the symmetry properties of the equilibrium distribution given as 
p(°°) ( ri i r2 ) — p(°°) ( r2) we can immediately see that the flux in this state vanishes: 
S(°°\R, R) = 0. This is of course a result expected for an asymptotic regime in the 
infinite medium. In order to describe states with a non-zero flux, which are most 
relevant to experimental situations, we have to consider the correlation function at a 
finite distance from the sources or in a finite-sized medium. In this case, we cannot 
expect that the same modes ^ will diagonalized the correlation function, but if the 
deviations from the equilibrium state are small we can assume that non-diagonal 
elements decay very fast when one moves away from the main diagonal of the density 
matrix. Then, instead of trying to find the true diagonalizing states we will follow 
the same approach as in the standard transport theory in statistically homogeneous 
medium and convert the non-diagonal density-matrix into a slowly changing in the 
real space function. The field-field correlator in this representation takes the form of 

p(r 1 ,r 2 )=J2 I kMR)^kMri)n^(r 2 ), (56) 

k,m 

where R = (t*i + r 2 )/2 and Ik,m(R) are smooth functions of the coordinate R. The 
functions lk,m(R) can be interpreted as specific intensities of the respective modes 
^k.m- In order to show it let us calculate the average energy density and the average 
Poynting vector. Using representation ([56]) in Eqs. ([32]) and (|3T[) we obtain 

(w(R)) = }^2k,m(R)Wk,m(R), 

(5?) 

(S(R)) = ^ Ik,m {R)Sk,fh (R) ) 
k.rh 

where we have used Eqs. (|35j) and (|40|) to define formally the energy density, Wk,m 
and the Poynting vector, Sh,m, corresponding to the states ^>k.m- Equations ([57)) 
agree with the intuitive concept of the specific intensity justifying our identification 
for functions Ik,m- 

In order to derive an equation for the specific intensities, we calculate the matrix 
element of both sides of ([24]) between the states ^>* k+q / 2 fi and ^k-q/ 2 ,a and then use 
the assumption of smooth Tk,m(R) as expressed by f|C8j) . Calculations based on the 
separation of length scales (see details in |Appendix CP lead to the radiative transfer 
equation in the form 

ijj Vm(k) • V B.Zk,m(R) = Im[Em i7fl (k,q)l q ^{R). (58) 

^ n 

Here and below the matrix elements for the self-energy and the irreducible vertex 
in the basis provided by the states are related to the matrix elements defined in 
Eqs. ([6]) and ([20]) through transformation (|46|) . 

In equation (|58|) we have introduced several important quantities. The velocity 
Vm(k) determines the direction of the "ray" corresponding to the mode ^>k,m and is 
expressed in terms of the group velocity of the photonic crystal modes as 

Vjn(k) = - V \B^ n (k)\ 2 LO n {k)v n (k). (59) 
0J 
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The spatial decay of the initial distribution Xk,m(R) concentrated in a single mode 
towards equilibrium is quantified by the mean-free-path 

^(fe) = - Im f"f fc f ] . (60) 

As follows from the optical theorem, the decay is a result of the redistribution of 
energy between the modes specified by different points in the first Brillouin zone and 
different band numbers. Making use of representation ((20)) the scattering kernel, which 
describes the redistribution, can be shown to have the form 

o-fn,n{k, q) = Np [ ^\k)U™-™(k 1 k, q, q). (61) 

Initially, in relatively small samples or near the boundary of a big sample, the 
distribution of energy can be arbitrary complex. For example, almost all energy can be 
concentrated in only few modes. The distribution changes along the sample according 
to ([58]) and eventually 2k,m{R) tends to p k °°l, which cancels both sides of ([58)) . This 

confirms our interpretation of p^°X as the asymptotic equilibrium distribution. 

If we can neglect the interband coupling and use the on-shell approximation, which 
is justifiable in the limit of weak and short-scale disorder, the picture significantly 
simplifies since in this approximation B rTl , n (k) become just identity matrices and the 
density matrix can be approximated by (|52[) . This means that the transport is provided 
by the modes of the ideal photonic crystal corresponding to the frequency ui and the 
specific intensity Xk,m has the transparent meaning of the intensity of these modes. 
The spatial distribution of the specific intensity is governed by the same equation ([55| 
with the velocity Vm(k) substituted by the respective group velocity Vm(k) — » v m (k) 
as follows from ([59| . 

In order to illustrate the specific features of the radiative transfer equation in 
disordered photonic crystals, we calculate the scattering kernel av^^g) f° r a 2D 
structure made of dielectric discs in the limit of weak disorder. We approximate the 
intensity propagator by the sum of the ladder diagrams and additionally assume the 
(5-correlated disorder when K(ri,r2) — V 2 S(ri — rs)5(ri — r2), where r$ denotes 
points lying in a thin shell near the surface of the ideal disc [see |Appendix A| . In this 
case the specific intensity is concentrated at the equifrequcncy surfaces F m (u>) so that 
it can be presented as Ik, m {R) = ik,m{R)5 [uj 2 — w^(fe)] . The amplitudes Ik, m {R) 
satisfy the radiative transfer equation of the same structure as (|58l) with the scattering 
amplitude between the states on the equifrequency surface given by 

<T m , n (k,q) = „ / dr |u fc . m (r)| 2 |Mq,„(r)| 2 , (62) 

Vn(q)V 2 J s 

where the integration is performed in a thin shell near the surface of the ideal disc. 
We show the dependence of a m . n (k, q) on the direction of the final Bloch wave vector 
in Figure ^Bp for two cases when the frequency ui is well below the frequency of the 
fundamental gap loq (log ~ 0.28 for the structure used in the numerical calculations) 
and when uj>uog- As one can see, for low frequencies the scattering amplitude is close 
to isotropic (the variation of the scattering amplitude is not visible in this scale), which 
is the consequence of slow spatial variation of the Bloch amplitudes Uk,m(j") and the 
shape of the equifrequency surface close to spherical. However, when u> approaches log 
both the spatial variation of the Bloch amplitudes and the shape of the equifrequency 
surface become nontrivial resulting in strongly anisotropic scattering. The interesting 
feature of this dependence is that the scattering rates between the states, which belong 
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Figure 2. (Colour online) (a) The band structure of the ideal photonic crystal 
with the same parameters as in Figure [T] Only two bands used in calculations for 
Figs. Q] [2p and [3] are shown. The horizontal lines correspond to the frequencies 
Lua/2n = 0.21 and 0.41. (b) The scattering amplitude cr(k, q)/w i V 2 between the 
states belonging to the first (uia/2ir = 0.21, the internal points) and the second 
(uia/2n = 0.41, the points corresponding to higher values) bands as a function of 
the angle between the vectors fc and q. The direction of the initial Bloch wave 
vector q is taken along the TX direction for both frequencies, thus the angle in 
(b) is counted from the YX direction. 



to the same band (the points corresponding to Loajlix = 0.41 near the horizontal line), 
are noticeably higher than the scattering amplitudes connecting the states in the 
different bands (points near the vertical axis). It should be noted that since the gap 
is not complete the scattering amplitude remains finite for all directions. 

The expression for the mean- free-path i m (k) in terms of the scattering amplitude 
o~m,n{k,q) is found from the requirement for p' 00 ' to solve the radiative transfer 
equation, 

£m( k ) = "7m / dgcr mi „(fc,<z), (63) 

where the integral runs over the respective equifrequency surfaces. Substituting (|62|) 
into this expression we find ^~ 1 (fe) in terms of the Bloch amplitudes. The same result 
can be obtained from (|60[) using for the self-energy the same approximation as for 
deriving ([62]). 

The evaluation of the mean-free-path for the same frequencies as in Figure [2] is 
presented in Figure El It is seen that for the relatively low frequencies the mean- free- 
path only weakly depends on the direction of the Bloch wave vector thus reproducing 
the well-known results for statistically homogeneous media. In a vicinity of the partial 
band gap the anisotropic scattering manifests in strong directional dependence of the 
mean- free-path. This is especially evident for an intermediate frequency ua/2n = 0.31, 
for which near the TM-direction the mean- free-path drops by the factor of 5. 
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Figure 3. (Colour online) The directional dependence of the mean-free-path 
£(u})uj 4 /V 2 for frequencies uja/2n = 0.21 (curve 1), 0.31 (curve 2) and 0.41 (curve 
3). The angle is counted from the TX-direction. 



3.4- The diffusion equation 

The right-hand side of (f55|) describing spatial variations of the specific intensities 
consists of the sum of two terms. In order to appreciate their physical significance let 
assume that at small distances from the source only single mode is excited so that the 
specific intensity can be presented in the following form 

l k , m oc l(R)S mtmo S(k - fc ). (64) 

The first of the right hand terms in (|58[l describes exponential spatial decay of this 
specific intensity due to scattering and redistribution of the energy among other modes. 
At small distances this term dominates and solution of (|58|) can be presented in 
the same form as (|64p with exponential spatial dependence. However, at distances 
from the source exceeding the mean-free-path this term diminishes and the second, 
integral term, starts determining the spatial distribution of the specific intensities, 
which is no longer exponential. One can describe this situation at asymptotically 
large distances taking into account that in this limit the density matrix /c(Ti,r 2 ) 
must be close to its limiting value p' 00 ) (ri, r-z) so that one can present p{r\,r2) 
using an appropriate expansion near pj^i- Drawing an analogy with the derivation 
of the diffusion approximation in the standard case [29l [42l [27] and assuming that 
the photonic crystal has the centre of symmetry we present specific intensities in the 
following form 

IhA R ) = —7TT [W-'w^R) + n m (fe)T- 1 5 d ( J R)] , (65) 

where Wd(R) and Sd(R) are assumed to be slowly changing functions of R, while W 
and T are scalar and tensor respectively, which do not depend on spatial coordinates. 
With an appropriate choice of these quantities one can show that Wd(R) and Sd(R) 
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in ([63)1 can be presented as 



w d (R) = — [ drw(r,r), 

1 f (66) 
S d (R) = - drS(r,r), 

17 -/V(H) 



where the integration is performed over the elementary cell with the coordinate 
R; w(ri,r 2 ), and S(r 1 ,r 2 ) are defined by Eqs. ([3"T]) and ([55]). respectively. These 
expressions clarify the physical meaning of Wd(R) and Sd(R) as long scale envelopes 
of the averaged energy density and the flux, respectively. W and T in ([65]) are found 
to be 

W=^X(u), T = f Akv m {k)®v m {k). (67) 

where ® denotes the tensor product, {v®u)ij — ViUj and A(w) is the density of states 
of the photonic crystal defined as A(oj) = ~ 

The equations with respect to u;<i(-R) and Sd(R) can be derived from the radiative 
transfer equation. Summation over all states on the respective equifrequency surface 
yields the energy conservation law 

V • S d (R) = 0. (68) 

Multiplying both sides of l[58p by v m (k) and summing over all states we obtain 

VDViu d (.R) = 0, (69) 

where the diffusion tensor D(uj) up to a constant factor is defined by D(w) = TM _1 T 
with the current relaxation kernel 1431 



Wdfe [C(*)-E / dg °^ ( ^ ) g m (fc)»g n (q) 



(70) 



It is important to note that despite the scalar character of the initial Helmholtz 
equation ([T]) the diffusion in photonic crystals is characterized by a diffusion tensor 
rather than by a single diffusion constant. In particular, if the periodic modulation 
is characterized by a rectangular (noncubic) elementary cell then the eigenvalues of 
the tensor corresponding to different principal directions will be different resulting in 
anisotropic diffusion. 

This diffusion anisotropy, however, appears only in low symmetry photonic 
crystals. If the point symmetry group of a particular structure is such that irreducible 
representations of the full rotation group remain irreducible for the point group of 
the crystal, the diffusion tensor is reduced to a scalar form. Indeed, crystals with 
such symmetries, e.g. crystals with square and hexagonal lattices in 2D and cubic 
and fee lattices in 3D, do not allow non-trivial tensors of the second rank [U] In this 
case the diffusion, despite of the strong directional dependence of the scattering kernel 
o~ mn (k, q), is isotropic and is characterized by the diffusion constant 

D= m_ %) (71) 

d\(u) 1 - (cos(0)> 1 ' 

where d is the dimensionality of the problem, F(u>) = J2 m J dk 1 is the total area of 
the equifrequency surface, 

^M = ^)E / dfc C(fc) (72) 



Radiative energy transfer in disordered photonic crystals 



20 




Figure 4. (Colour online) The dependence of the effective mean-free-path 
£(u))/V 2 (dashed line, squares) on frequency in log- log scale. The solid line 
(circles) shows the dependence oc ui 3 (this asymptotics for low frequencies is due 
to the dimensionality of the problem). 



and 

<«*(*)> ^ #4 E / dfed ? ^ ■ ^*>- ( ?3 ) 

The expression for the diffusion coefficient has the same general structure as for 
statistically homogeneous media with the transport velocity Ve(uj) — F{lS)/\(u>) 
and the "transport mean-free-path" £(uj)/(1 — (cost?)). The transport velocity is 
independent on the details of the disorder distribution owing to the ^-functional 
approximation for the correlation function K(ri,r2), which is valid if the frequencies 
corresponding to morphological resonances at the typical inhomogcncity size is much 
higher than the range of frequencies under consideration [J5J USJ H7] . Figure [5] shows 
the frequency dependence of £(u>) for a disordered photonic crystal with the square 
lattice. It should be noted that the minimum of £(u>) is reached at the band edge 
in the TX-direction (see Figure [2^) and is due to the significant drop of the group 
velocity when the equifrequency surface touches the boundary of the first Brillouin 
zone. 

It should be emphasized, however, that the possibility to introduce a single scalar 
diffusion constant and, respectively, to describe transport by a single transport mean- 
free-path is not a a result of the disorder destroying the effect of periodicity but rather 
is the consequence of the underlying symmetry of the photonic crystal. In the case of 
low symmetry structures the diffusion is described by a tensor, and is characterized 
by different effective mean-free-paths for different principle directions of the diffusion 
tensor. 

4. Conclusion 

In this paper we developed a systematic approach based on multiple scattering analysis 
to the theoretical description of incoherent transport properties of disordered photonic 
crystals in the steady state regime. The main difficulty in developing such a theory 
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using a standard plane-wave based multiple scattering approach lies in the separation 
between waves scattering due to periodic modulations of the dielectric constant 
and scattering due to random deviations from the periodicity. One might think 
that this difficulty is readily resolved by incorporating modes of the ideal photonic 
crystal instead of plane waves into the standard transport theory. However, such a 
straightforward approach fails because of the inherent presence of two spatial scales, 
the period of the crystal and the main-free-path, both of which has to be retained 
in the theory. As a result, a formally constructed Wigner distribution, which is the 
main technical tool in deriving the radiative transfer equation in the regular case of 
statistically homogeneous medium, loses its smooth spatial dependence and, hence, 
its methodological usefulness. 

We showed that these difficulties can be overcome if one incorporates the 
concept of the radiative transfer as a incoherent process into the foundation of the 
theory. In order to achieve this, we showed that the field-field correlation function, 
p('Ti, r 2) — (E(ri)E*(r2)), has all the properties of a density matrix defined in 
quantum statistics, and hence, can be interpreted as such. This interpretation allowed 
us to separate coherent and incoherent contributions to the transport and eventually 
obtain radiative transfer and diffusion equations for media with periodically modulated 
average dielectric function. 

We found an exact asymptotic solution, p(°°>, of the Bethe-Salpeter equation (Ti2"|) , 
which describes the limiting form of the correlation function in an infinite medium 
asymptotically far away from the sources. This function determines the asymptotic 
distribution of the intensity of the field, which was shown to be highly spatially 
non-uniform and anisotropic. This result should be contrasted with the case of 
statistically homogeneous medium, in which the asymptotic distribution of intensity 
is spatially uniform. It is important that the properties of the intensity distribution 
even inside an infinite photonic crystal are determined by the normal modes of the 
underlying periodic structure. This result casts serious doubts on the assumption 
made in fiefs. [HI [16] that the field distributions deep inside a photonic crystal and a 
statistically homogeneous media do not differ from each other, and that the anisotropy 
of the light emerging from photonic crystals is formed within a narrow layer at the 
boundary. 

The asymptotic character of p^ 00 -* is expressed by the absence of flux in this 
state. Using analogy with the statistical physics we can interpret p(°°' as an 
equilibrium distribution. In order to describe actual energy transfer, one needs 
to consider small deviations from equilibrium. We derived generalized forms of 
radiative transfer and diffusion equations valid for disordered photonic crystals and 
obtained general expressions for the scattering cross-section, mean free path, and 
the diffusion coefficient. As an example, we calculated the cross-section for one 
particular model of a photonic crystal and demonstrated how the underlying periodic 
structure effects disorder induced scattering of photonic modes. In particular, we 
found that the scattering cross-section describing the redistribution of the energy 
between modes becomes highly anisotropic at high frequencies. The mean-free-path 
of a particular mode, which describes the rate of redistributing the correspondent 
specific intensity, also depends strongly on the direction of the respective Bloch vector. 
The numerical calculations show that for a chosen frequency near the (partial) band 
gap the mean-free-path may vary almost by the order of magnitude. Such significant 
variation presents especial interest from the perspective of the problem of the Anderson 
localization. 



Radiative energy transfer in disordered photonic crystals 



22 



The interesting feature of the transport in the disordered photonic crystals is 
that even if scattering between different modes may be highly anisotropic it does not 
necessarily imply an anisotropic diffusion. Indeed, if the point symmetry of the crystal 
is sufficiently rich, e.g. in 2D this is rotational symmetry of the third order or higher, 
then all tensors of the second rank describing the intrinsic macroscopic properties 
are proportional to the unit tensor. That is the diffusion is necessarily isotropic and 
is characterized by a single diffusion constant. The mean-free-path appears in this 
constant averaged over the respective equifrequency surface. 

The derivation of the transport equations can be generalized to describe 
the radiative transport in more general situations than considered in the paper. 
Principally, the equilibrium distribution given by p(°°' holds in the presence of 
the boundary as well. The latter is accounted through the boundary conditions 
determining Green's functions. However, whether the respective density matrix can 
be described near the boundary as a long-scale perturbation of or not presents a 
special interesting problem. 

Finally, we would like to comment on the possibility to extend the consideration 
provided in the paper to a non-steady-state regime. By a direct analogy with the 
standard case, a slow dynamics in the diffusion regime can be accounted for by 
assuming Wd in (|65p to be time dependent. However, a rigourous theory, which 
would yield this limit is yet to be developed. From the perspectives of the presented 
consideration, the main formal obstacle is the fact that Pu lt u? is generally not non- 
negatively defined if u>i ^ u>2- This problem, of course, exists and is recognized in 
the standard theory of multiple scattering in homogeneous media as well. In order 
to ensure positivity of the specific intensities, it was suggested to consider a coarse- 
grained Wigner distribution [48j . A possibility to generalize our approach to time- 
dependent correlation functions and their possible relation with the density matrix 
formalism is an open problem. 
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Appendix A. A model of Gaussian inhomogeneities in photonic crystals 

A Gaussian random field Ae(r) can be defined in two equivalent ways. The functional 
definition is based on the requirement of an integral 



to be a Gaussian random variable for an arbitrary function f(r). The statistical 
definition states that a multi-point correlation function of a zero-mean Gaussian field 
is expressed in terms of the two-point correlation function 



where the summation is taken over all possible pairings of the points r% , . . . , r2 n and 




(A.1) 



(Ae(n) . . . Ae(r 2n )) = £ K(r a , r fJ )... K{ 



(A.2) 



pairs 



K(n,r 2 ) = <Ae(n)Ae(r 2 )>. 
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Property of integral (|A. 1|) being a Gaussian random variable is consistent with the 
intuitive idea about a Gaussian random field, while property (|A.2[) is convenient from 
the technical point of view for developing the perturbation theory [55]. Both these 
properties, however, might seem somewhat artificial from the perspective of disordered 
photonic crystal. Here we would like to present a simple model of the disorder in 
a photonic crystal that naturally results in simulation of the inhomogeneities by a 
Gaussian random field. 

Let the ideal structure be constituted by spheres with the dielectric constant e\ at 
the sites of a periodic lattice. Furthermore, we assume that the disorder in this system 
is constituted by slight variations of the size of the spheres, while their positions are 
fixed. Thus, we can represent the spatial modulation of the dielectric function in the 
disordered crystal in the form 

e(r) = e +J2^R(r + R), (A.3) 

R 

where R are the lattice vectors, eo is the background dielectric constant, and 
StRir + R) is the deviation from the background value in the elementary cell with the 
coordinate R. By the construction <5ejt(r) = t\ — eo inside the sphere and it equals 
to elsewhere. We represent the spatial distribution of the dielectric function in the 
form adapted in ((T|), i.e. e(r) — e(r) + Ae(r), where e(r) describes the distribution in 
the ideal structure and 

Ae(r) = O^O + R ) ~ 5e o(r + R)) (A.4) 

R 

with Seo(r) being the deviation from the background value in the ideal structure. 

Clearly Ae(r) is not zero in shells near the boundaries of the spheres constituting 
the ideal structure. It is positive in the spheres with the radius larger than the radius 
of the ideal spheres and is negative in smaller spheres. We simulate this function by 
presenting it in the form 

Ae(r) =^Ae H w(r- + i?), (A.5) 

R 

where u(r) is a non-random function different from zero in a thin shell near the 
boundary of the ideal sphere and the random amplitudes Aejj have the Gaussian 
distribution with zero mean value and are independent in different elementary cells. 
Using (|A.5|) in the integral in (jA.lj) one has a Gaussian random variable and, hence, 
Ae(r) is the Gaussian random field. 

This model, obviously, can be generalized by allowing the random variables Aej? 
to be a (usual) Gaussian random field rather than a constant inside each elementary 
cell. Such model of the disorder would account not only for the size dispersion of the 
spheres but also for the roughness of their surfaces. Assuming that Ae r in different 
elementary cells are not correlated we obtain 

K(n,r 2 ) = (Ae Pl Ae Pa >u(ri)«(r 2 ) (A.6) 

if ri and r-i are situated inside the same elementary cell and K(r\,T2) = otherwise. 
For calculations in the main text we have used the simplest model of spheres with 
uniform sizes and with centers at the sites of an ideal lattice but with rough surfaces. 
Thus we take u(r) = 1 in a thin shell around the ideal sphere (respectively, the circle 
in 2D) and 

(Ae ri Ae r2 ) = V 2 5(n-r 2 ). (A.7) 
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Appendix B. Symmetries of the averaged Green's function, the 
self-energy and the irreducible vertex 

Leaving the problem of convergence of the perturbational series aside, the symmetry 
properties of the averaged Green's function (G(ri,r 2 )), the self-energy £(ri,r 2 ), and 
the irreducible vertex U(ri, r 2 ; r[, r 2 ) can be proven on the diagram by diagram basis 
[26] , We demonstrate the typical line of arguments by proving that if the correlation 
function of inhomogeneities is invariant with respect to lattice translations then so is 
the averaged Green's function. 

Any diagram in the perturbational expansion of (G(ri,r 2 )) has the form of a line 
with 2n internal (n > 0) and 2 terminating vertices. The internal vertices divide the 
line on 2n + 1 segments, which are the graphical representations of the non-perturbed 
Green's function of the ideal structure. The internal vertices are connected pair-wise 
by the lines of the correlation function. 

Clearly the value of a diagram with fixed values of the coordinates corresponding 
to the vertices does not change if we shift the coordinates of all points (internal and 
terminating) by the vector of the lattice translation. Indeed, the Green's function 
lines and the lines corresponding to the correlation functions do not change because 
of the translational invariancc of these functions. Furthermore, we note that in order 
to obtain the contribution of the diagram into (G(ri,T"2)) we need to integrate with 
respect to coordinates of the internal vertices. Finally, we observe that the uniform 
shift of the coordinates of all vertices produces the diagram, which corresponds to the 
perturbational series for (G(ri + R, r 2 + R)). 

This line of arguments proves the translational invariance of the averaged Green's 
function, the self-energy and the irreducible vertex. In order to prove the reciprocity 
of the Green's function and the self-energy 

(G(n, r 2 )> - (G(r 2 , n)) , £(n, r 2 ) = £(r 2 , n), (B.l) 

we need the version of the arguments sketched above, which includes the directionality 
of the lines constituting the diagram. The reciprocity would follow then from 
invariance of the diagram with respect to the reversion of the directionality of all 
lines. For the irreducible vertex we additionally need to take into account that the 
correlation function connecting upper and lower lines is mapped to itself upon the 
reverting. Thus, we have 

U(n, r 2 ; r[,r' 2 ) = U(r[,r' 2 ; r u r 2 ). (B.2) 
We complete the consideration of the symmetry properties by proving the 
invariance of the self-energy with respect to the point symmetries of the photonic 
crystal. More specifically, we show that if the correlation function transforms according 
to the identity representation of the respective group, then so does the self-energy. 
The proof is based on the property of the unperturbed Green's function Go(t*i, r 2 ) to 
transform according to the identity representation 

TG (ri,r 2 ) EGofT^riJ-'rj) = G (n, r 2 ), (B.3) 

where T is an element of the point symmetry group. The proof of the invariance of 
the self-energy uses the same arguments as above with the only difference that now 
instead of shift or inversion operators we act with T on all points of the diagram. 

It follows from this consideration that the eigenfunctions of the Dyson equation 
can be classified according to the irreducible representations of the group of the point 
symmetries of the ideal structure. We use this fact in the paper when we discuss the 
effect of the disorder on the degenerate points in the spectrum of the photonic crystal. 
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Appendix C. Separation of spatial scales in close-to-equilibrium regimes 

The main assumption used for the derivation of the radiative transfer equation and 
the diffusion equation in the paper is that the spatial scales of variation of the specific 
intensity or the envelope of the energy density and of the functions describing the 
modes of the ideal photonic crystal can be separated. Physically this assumption 
implies that the slowly varying functions do not lead to coupling between different 
Bloch modes. As a result, one can neglect the smooth functions while calculating the 
respective scalar products. 

In order to give a formal expression for the separation of length scales we consider 
a photonic crystal mode ^k.mir) modulated by a smooth function f(r) 

h k , m {r) = /(r)* fe , m (r). (C.l) 

The smoothness of the function f(r) can be quantified by (p)j and {p 2 ) p where 

(p) f = J dppf(p) (C.2) 

and f(p) is the Fourier image of f(r). The limit of smooth f(r) can, thereby, be 
formalized as (p) j — * and (p 2 ) . — » 0. In this limit the weighted scalar product of 
hk,m(r) with the mode ^g,n(r) is written as 



(*q,n, h h , m ) = J dr e{r)% in {r)h h>m {r) = f(k - q)U n , m (q, fe), (C.3) 



where 



U n . m (q,k) = / dre(r)u* (r)u k , rn (r). (C.4) 
JV 

The inverse Fourier transform of (|C.3|) with respect to k — q gives for m = n 

J dp (* fe _ p/2 , m , h k+p/2 , m ) e-^ r = f(r) + 0(/'(r)), (C.5) 

where the remainder is small together with the derivatives of /(r). 

Similarly the notion of the smooth modulation can be applied for functions of 
the form 

p(ri,r 2 )=^T fe , m [(r 1 + r 2 )/2]* fe>ro (r 1 )^ m (r 2 ), (C.6) 

k,m 

where Ik,m{ r ) are smooth functions. The weighted matrix element of p(ri,r2) 
between * p+ q/ 2 ,„ and *p_g/ 2 , n is 

(^p+q/2,n\p\^p-q/2,n) =^2^p,m(q)Um,n(P + q/2,P)Un,m(P,P ~ <?/2). (C.7) 

m 

The inverse Fourier transform with respect to q gives 

J dqe- iq R (*p+,/2,n| P |*p-«/a,n) = l P ,m(R) (C8) 

up to terms vanishing with the derivatives of Ik, m (r). 

Additionally, for derivation of the radiative transfer equation one needs to 
compare the scales related to the irreducible vertex. Analyzing the diagrammatic 
expansion for U(r 1 ,r2;r' 1 ,r' 2 ) one can see that there are two typical spatial scales 
determining the decay of the irreducible vertex with the distance between the first and 
the second pairs of the points. The "short" scale is proportional to the correlation 
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radius of the inhomogeneities. This spatial decay is characteristic for ladder and 
maximal crossed diagrams approximations. Assuming that the correlations vanish at 
the scale of the elementary cell one can approximate 



This approximation has been used for derivation of (|58p . 

There are additional terms in the diagrammatic expansion of the irreducible 
vertex that go beyond the ladder and maximally crossed diagrams approximations 
and lead to the spatial decay on the scale of the mean-free path. These terms would 
result in a non-local term in the radiative transfer equation. The effect of the non- 
local scattering on the radiative transfer in disordered photonic crystals will be studied 
elsewhere. 
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